#include "utils.h"

#include <stdlib.h>

#include <math.h>
#include <sys/time.h>
#include <limits.h>

double standard_deviation(double *d, int n, double mean)
{
    int i;
    double x = 0.0;

    for( i=0; i< n; ++i )
        x += (d[i]-mean)*(d[i]-mean);

    return sqrt(x/n);
}

const unsigned long A = 1664525, C = 1013904223;
inline unsigned long random_ul(unsigned long x){
    return x*A+C;
}

inline double random0_1(unsigned long *x){
    return (*x = random_ul(*x))/(double)ULONG_MAX;
}

void initBoxMullerState(struct BoxMullerState* state)
{
    state->random.__init = 0;
    state->useLast = 0;

    struct timeval now;
	gettimeofday(&now, NULL);
	state->random.__x[0] = now.tv_usec;
}

double boxMullerRandom(struct BoxMullerState* state)
{
    if (state->useLast)
    {
        state->y1 = state->y2;
        state->useLast = 0;
    }
    else
    {
        do
        {
            drand48_r(&state->random, &state->x1);
            state->x1 = 2.0 * state->x1 - 1.0;
            drand48_r(&state->random, &state->x2);
            state->x2 = 2.0 * state->x2 - 1.0;
            state->w = state->x1 * state->x1 + state->x2 * state->x2;
        }
        while (state->w >= 1.0);

        state->w = sqrt((-2.0 * log(state->w)) / state->w);
        state->y1 = state->x1 * state->w;
        state->y2 = state->x2 * state->w;
        state->useLast = 1;
    }

    return state->y1;
}
